## Evaluation of Wahlomat experiment -----

# load packages
source("packages.r")


## data import ------------------------------

load("../data/yougov_panel_ger_wide_wahlomat_prep.RData")



### Model main effects -----------------------------
# setup instrument, treatment, covariates
assignment <- "wahlomat_encouragement"
treatment <- "wahlomat_use_combined"
covar_names <- c("gender.1", "age.1", "educ.1", "hhincome.1", "polinterest.1", "leftright.1")
#covar_names <- c("gender.1", "age.1", "educ.1", "hhincome.1", "polinterest.1", "leftright.1", "voted_2013_yesno", "civic_knowledge_2")


## Generic function to identify effects under naive, ITT, CACE model ----

naive_itt_cace_estimation <- function(dv, label, scale, iv_robust=TRUE, model_out = "none") { 
  # naive model based on observed treatment (disregarding assignment)
  observed_formula <- paste(dv, " ~ ", paste(c(covar_names, treatment), collapse=" + "))
  model_observed <- lm(observed_formula,  data = waves_df_wide)
  summary(model_observed)
  
  # ITT model based on encouragement assignment
  itt_formula <- paste(dv, " ~ ", paste(c(covar_names, assignment), collapse=" + "))
  model_itt <- lm(itt_formula,  data = waves_df_wide)
  summary(model_itt)
  
  # IV model to identify CACE
  iv_formula <- paste(dv, " ~ ", 
                      paste(c(covar_names, treatment), collapse=" + "), "|", 
                      paste(c(covar_names, assignment), collapse=" + "), sep = "") %>% as.formula
  if (iv_robust==FALSE){ model_cace <- ivreg(iv_formula,  data = waves_df_wide) }
  if (iv_robust==TRUE){ model_cace <- iv_robust(iv_formula,  data = waves_df_wide) }
  summary(model_cace, diagnostics = TRUE)
  
  # computing F-stat from first stage
  restricted <- lm(paste(treatment, "~", 
                      paste(covar_names, collapse=" + "),
                      sep=""), data=waves_df_wide)
  unrestricted <- lm(paste(treatment, "~", 
                       paste(c(covar_names, assignment), collapse=" + "),
                       sep=""), data=waves_df_wide)
  fstat <- anova(restricted, unrestricted)$F[2]
  message("F-test for weak instrument in 1st stage = ", round(fstat,2))

  # get estimates
  model_observed_predict <- broom::tidy(model_observed)
  model_itt_predict <- broom::tidy(model_itt)
  model_cace_predict <- broom::tidy(model_cace)
  
  model_effects <- rbind(filter(model_observed_predict, str_detect(term, "wahlomat")),
                         filter(model_itt_predict, str_detect(term, "wahlomat")),
                         filter(model_cace_predict, str_detect(term, "wahlomat"))[
                           names(model_itt_predict)
                         ])
  
  model_effects$model <- c("Observed", "ITT", "CACE")
  model_effects$dv <- label
  model_effects$scale <- scale
  
  model_effects <- mutate(model_effects, 
                          lo95 = estimate - 1.96*std.error,
                          hi95 = estimate + 1.96*std.error,
                          lo83 = estimate - 1.4*std.error,
                          hi83 = estimate + 1.4*std.error)
  if(model_out == "none") {
    model_effects
  }else if(model_out == "CACE"){
    model_cace
  }else if(model_out == "ITT"){
    model_itt
  }else if(model_out == "Observed"){
    model_observed
  }else if(model_out == "First stage"){
    unrestricted
  }
}


## generic function to run texreg
texreg_models <- function(){
  texreg(list(
    naive_itt_cace_estimation(dv, label, scale, model_out = "Observed"),
    naive_itt_cace_estimation(dv, label, scale, model_out = "ITT"),
    extract.iv_robust(naive_itt_cace_estimation(dv, label, scale, model_out = "CACE"), include.ci = FALSE, include.adjrs = FALSE)),
    digits = 3,
    custom.model.names = c("Observed", "ITT", "CACE"),
    custom.coef.names = c("Intercept", "Gender", "Age", "Education", "Household Income", "Political Interest", "Left-Right Ideology", "Wahlomat Use", "Wahlomat Encouragement"),
    caption = paste0("Effect of VAA usage/encouragement on ", tolower(label), "."),
    caption.above = TRUE,
    fontsize = "normalsize",
    include.adjrs = FALSE,
    single.row = TRUE,
    booktabs = TRUE,
    dcolumn = TRUE,
    use.packages = FALSE,
    label = paste0("tab:reg", dv), 
    file = paste0("../output/regtable-", dv, ".tex")
  )
}


#### 0) FIRST STAGE MODEL ---------------

dv <- "secondvote_3to4"
label <- "First stage"
scale <- "[0-1]"
model_first_stage <- naive_itt_cace_estimation(dv, label, scale, model_out = "First stage")
texreg(model_first_stage,
       digits = 3,
  custom.model.names = c("Wahlomat Use"),
  custom.coef.names = c("(Intercept)", "Gender", "Age", "Education", "Household Income", "Political Interest", "Left-Right Ideology", "Wahlomat Encouragement"),
  caption = paste0("First-stage effects of VAA encouragement and covariates on VAA usage"),
  caption.above = TRUE,
  fontsize = "normalsize",
  include.adjrs = FALSE,
  single.row = TRUE,
  booktabs = TRUE,
  dcolumn = TRUE,
  use.packages = FALSE,
  label = "tab:regfirststage", 
  file = paste0("../output/regtable-firststage.tex")
)



#### 1) MODEL MAIN RESULTS ---------------

## Change in vote choice from wave 3 to 4 -----------------------

waves_df_wide$secondvote_3to4 <- waves_df_wide$secondvote.3 != waves_df_wide$secondvote.4
dv <- "secondvote_3to4"
label <- "Switched vote choice"
scale <- "[0-1]"
table(waves_df_wide[,dv])
prop.table(table(waves_df_wide[,dv]))
model_effects_follow_vote_switch <- naive_itt_cace_estimation(dv, label, scale)
power.prop.test(n=500, p1=0.79, power=0.80, sig.level=0.05)
texreg_models()

model_out <- naive_itt_cace_estimation(dv, label, scale, model_out = "CACE", iv_robust = FALSE)
model_out$nobs # number of observations
tabyl(model_out$model$wahlomat_use_combined) # number of (validated) VAA users


## Vote/turnout (reported) -----------------------

dv <- "voted_2017_yesno"
label <- "Turnout"
scale <- "[0-1]"
table(waves_df_wide[,dv])
prop.table(table(waves_df_wide[,dv]))
model_effects_voted_2017_yesno <- naive_itt_cace_estimation(dv, label, scale)
power.prop.test(n=500, p1=0.90, power=0.80, sig.level=0.05)
texreg_models()

model_out <- naive_itt_cace_estimation(dv, label, scale, model_out = "CACE", iv_robust = FALSE)
model_out$nobs # number of observations
tabyl(model_out$model$wahlomat_use_combined) # number of (validated) VAA users



## Knowledge on party positions -----------------------

dv <- "wahlomat_perceivedpos_correct_avg"
label <- "Issue knowledge" # "Knowledge on parties'\nissue positions"
hist(waves_df_wide[,dv])
scale <- "[0-1]"
model_effects_wahlomat_perceivedpos_correct_avg <- naive_itt_cace_estimation(dv, label, scale)
texreg_models()

model_out <- naive_itt_cace_estimation(dv, label, scale, model_out = "CACE", iv_robust = FALSE)
model_out$nobs # number of observations
tabyl(model_out$model$wahlomat_use_combined) # number of (validated) VAA users



# interpretation: increase in % correct from 62.3% to 66.2%
aggregate(wahlomat_perceivedpos_correct_avg~wahlomat_use, 
          waves_df_wide, FUN=mean)
waves_df_wide$recoded_interest.1 <- factor(ifelse(waves_df_wide$polinterest.1 %in% c(1,2), "Low",
                                           ifelse(waves_df_wide$polinterest.1 %in% c(4,5), "High",
                                                  "Moderate")),
                                           levels=c("Low", "Moderate", "High"))

aggregate(wahlomat_perceivedpos_correct_avg~recoded_interest.1, 
          waves_df_wide, FUN=mean)
# effect size is roughly the same size as the gap between moderate interest and high interest
# and 50% of the gap between low and high interest


## Change in news consumption -----------------------
dv <- "log(news_count+1)"
label <- "News consumption (news\nURLs visited)"
scale <- "[log count]"
hist(log(waves_df_wide$news_count+1))

model_effects_posts_news <- naive_itt_cace_estimation(dv, label, scale) 
texreg_models()

model_out <- naive_itt_cace_estimation(dv, label, scale, model_out = "CACE", iv_robust = FALSE)
model_out$nobs # number of observations
tabyl(model_out$model$wahlomat_use_combined) # number of (validated) VAA users



## Facebook political use -----------------------
dv <- "I(political_posts_count>0)"
prop.table(table(waves_df_wide$political_posts_count>0))
label <- "1+ political posts\non Facebook"
scale <- "[0-1]"
model_effects_posts_facebook <- naive_itt_cace_estimation(dv, label, scale) 
texreg_models()


# interpretation: increase in % correct from 7.8% to 11.6%
aggregate(political_posts_count>0~wahlomat_use, 
          waves_df_wide, FUN=mean)
aggregate(political_posts_count>0~recoded_interest.1, 
          waves_df_wide, FUN=mean)
# effect size is roughly the same size as the gap between moderate interest and high interest
# and 50% of the gap between low and high interest

## Twitter political use -----------------------
dv <- "I(political_tweets>0)"
prop.table(table(waves_df_wide$political_tweets>0))
label <- "1+ political tweets\non Twitter"
scale <- "[0-1]"
model_effects_posts_twitter <- naive_itt_cace_estimation(dv, label, scale) 
texreg_models()


# interpretation: increase in % correct from 16.6% to 32%
aggregate(political_tweets>0~wahlomat_use, 
          waves_df_wide, FUN=mean)
aggregate(political_tweets>0~recoded_interest.1, 
          waves_df_wide, FUN=mean)
# effect size is roughly the same size as the gap between moderate interest and high interest
# and 50% of the gap between low and high interest



## Combine 1st set of model results: turnout, vote choice, knowledge  -------

model_effects <- rbind(model_effects_voted_2017_yesno,
                       model_effects_follow_vote_switch,
                       model_effects_wahlomat_perceivedpos_correct_avg
)

dat <- filter(model_effects, model == "CACE")
labels = paste(as.character(dat$dv), dat$scale)

model_effects <- transform(model_effects, 
                           model = factor(model,levels = c("Observed", "ITT", "CACE")),
                           dv = factor(dv, levels = unique(model_effects$dv), labels = labels))

p <- ggplot(model_effects, aes(model, estimate)) + 
  geom_pointrange(aes(ymin = lo95, ymax = hi95), shape = 19, fatten = 1, size = 1) + 
  geom_pointrange(aes(ymin = lo83, ymax = hi83), shape = 19, fatten = 1, size = 2) + 
  geom_hline(yintercept = 0) + 
  facet_wrap( ~ dv, scales = "free", nrow = 1) +
  theme_ipsum_rc() + 
  theme(panel.background = element_blank(), 
        axis.title.y = element_blank(), 
        axis.title.x = element_blank(), 
        axis.text = element_text(size = 12), 
        plot.caption = element_text(size = 8), 
        strip.text.x = element_text(size = 12, colour = "darkblue", face = "bold"),
        plot.margin=unit(c(0.1,0.1,0.1,0.1), "cm"),
        panel.spacing = unit(0.25, "cm"))
p
# save plot
ggsave(p, file = "../output/wahlomat-main-effects.png",
       width = 3.5, height = 1.25, units = "cm", dpi = 300, scale = 5)



## Combine 2nd set of model results: news consumption, facebook, twitter posts  -------

model_effects <- rbind(model_effects_posts_news,
                       model_effects_posts_facebook,
                       model_effects_posts_twitter
)

dat <- filter(model_effects, model == "CACE")
labels = paste(as.character(dat$dv), dat$scale)

model_effects <- transform(model_effects, 
                           model = factor(model,levels = c("Observed", "ITT", "CACE")),
                           dv = factor(dv, levels = unique(model_effects$dv), labels = labels))

p <- ggplot(model_effects, aes(model, estimate)) + 
  geom_pointrange(aes(ymin = lo95, ymax = hi95), shape = 19, fatten = 1, size = 1) + 
  geom_pointrange(aes(ymin = lo83, ymax = hi83), shape = 19, fatten = 1, size = 2) + 
  geom_hline(yintercept = 0) + 
  facet_wrap( ~ dv, scales = "free", nrow = 1) +
  theme_ipsum_rc() + 
  theme(panel.background = element_blank(), 
        axis.title.y = element_blank(), 
        axis.title.x = element_blank(), 
        axis.text = element_text(size = 12), 
        plot.caption = element_text(size = 8), 
        strip.text.x = element_text(size = 12, colour = "darkblue", face = "bold"),
        plot.margin=unit(c(0.1,0.1,0.1,0.1), "cm"),
        panel.spacing = unit(0.25, "cm"))
p
# save plot
ggsave(p, file = "../output/wahlomat-effects-news-social.png",
       width = 3.5, height = 1.25, units = "cm", dpi = 300, scale = 5)



#### 2) MODEL ADDITIONAL OUTCOMES (APPENDIX) ---------------

## Change in likelihood to vote -----------------------

dv <- "likelihoodvote_3to4"
label <- "Change in likelihood\nto vote"
scale <- "[1-10]"
table(waves_df_wide$likelihoodvote.3)
table(waves_df_wide[,dv])
model_effects_likelihoodvote_3to4 <- naive_itt_cace_estimation(dv, label, scale)

## Change in vote certainty -----------------------

dv <- "votecertainty_3to4"
label <- "Change in vote\ncertainty"
scale <- "[1-10]"
table(waves_df_wide$votecertainty.3)
table(waves_df_wide[,dv])
prop.table(table(waves_df_wide[,dv]))
model_effects_votecertainty_d <- naive_itt_cace_estimation(dv, label, scale)

## Mobilized to vote (2013-2017) -----------------------

dv <- "mobilized_2017"
label <- "Mobilized to vote"
scale <- "[0-1]"
table(waves_df_wide[,dv])
prop.table(table(waves_df_wide[,dv]))
model_effects_mobilized <- naive_itt_cace_estimation(dv, label, scale)

## Change in sympathy towards CDU ----------------------------

dv <- "scaloparties_cdu_3to4"
label <- "Change in sympathy\ntowards CDU/CSU"
scale <- "[1-11]"
table(waves_df_wide[,dv])
table(waves_df_wide$scaloparties_cdu.1)
model_effects_wahlomat_scaloparties_change_CDU <- naive_itt_cace_estimation(dv, label, scale)

## Change in sympathy towards SPD ----------------------------

dv <- "scaloparties_spd_3to4"
label <- "Change in sympathy\ntowards SPD"
scale <- "[1-11]"
table(waves_df_wide[,dv])
table(waves_df_wide$scaloparties_spd.1)
model_effects_wahlomat_scaloparties_change_SPD <- naive_itt_cace_estimation(dv, label, scale)

## Change in sympathy towards FDP ----------------------------

dv <- "scaloparties_fdp_3to4"
label <- "Change in sympathy\ntowards FDP"
scale <- "[1-11]"
table(waves_df_wide[,dv])
table(waves_df_wide$scaloparties_fdp.1)
model_effects_wahlomat_scaloparties_change_FDP <- naive_itt_cace_estimation(dv, label, scale)

## Change in sympathy towards Greens ----------------------------

dv <- "scaloparties_gru_3to4"
label <- "Change in sympathy\ntowards Greens"
scale <- "[1-11]"
table(waves_df_wide[,dv])
table(waves_df_wide$scaloparties_gru.1)
model_effects_wahlomat_scaloparties_change_GRU <- naive_itt_cace_estimation(dv, label, scale)

## Change in sympathy towards Left ----------------------------

dv <- "scaloparties_lin_3to4"
label <- "Change in sympathy\ntowards Left"
scale <- "[1-11]"
table(waves_df_wide[,dv])
table(waves_df_wide$scaloparties_lin.1)
model_effects_wahlomat_scaloparties_change_LIN <- naive_itt_cace_estimation(dv, label, scale)

## Change in sympathy towards AfD ----------------------------

dv <- "scaloparties_afd_3to4"
label <- "Change in sympathy\ntowards AfD"
scale <- "[1-11]"
table(waves_df_wide[,dv])
table(waves_df_wide$scaloparties_afd.1)
model_effects_wahlomat_scaloparties_change_AFD <- naive_itt_cace_estimation(dv, label, scale)


## Political interest ----------------------------

dv <- "polinterest_1to4"
label <- "Change in political\ninterest"
scale <- "[1-4]"
table(waves_df_wide[,dv])
model_effects_wahlomat_change_interest <- naive_itt_cace_estimation(dv, label, scale) 

## Political efficacy ----------------------------

dv <- "efficacy"
label <- "Political efficacy"
scale <- "[sd = 1]"
summary(waves_df_wide[,dv])
model_effects_wahlomat_efficacy <- naive_itt_cace_estimation(dv, label, scale) 

## Change in civic knowledge ------------

dv <- "civic_knowledge_2to5"
label <- "Change in\ncivic knowledge"
hist(waves_df_wide[,dv])
scale <- "[0-2]"
model_effects_wahlomat_civic_knowledge <- naive_itt_cace_estimation(dv, label, scale)

## Change in candidate recognition ---------------

dv <- "candidate_knowledge_2to5"
label <- "Change in\ncandidate recognition"
hist(waves_df_wide[,dv])
scale <- "[0-8]"
model_effects_wahlomat_candidate_knowledge <- naive_itt_cace_estimation(dv, label, scale)


## Political event items ---------------

dv <- "event_knowledge.5"
label <- "Event knowledge"
hist(waves_df_wide[,dv])
scale <- "[0-10]"
model_effects_wahlomat_event_knowledge <- naive_itt_cace_estimation(dv, label, scale)


## Combine model results; plot -------

model_effects <- rbind(model_effects_likelihoodvote_3to4,
                       model_effects_votecertainty_d,
                       model_effects_mobilized,
                       model_effects_wahlomat_change_interest,
                       model_effects_wahlomat_efficacy,
                       model_effects_wahlomat_scaloparties_change_CDU,
                       model_effects_wahlomat_scaloparties_change_SPD,
                       model_effects_wahlomat_scaloparties_change_FDP,
                       model_effects_wahlomat_scaloparties_change_GRU,
                       model_effects_wahlomat_scaloparties_change_LIN,
                       model_effects_wahlomat_scaloparties_change_AFD
                       )

dat <- filter(model_effects, model == "CACE")
labels = paste(as.character(dat$dv), dat$scale)

model_effects <- transform(model_effects, 
                           model = factor(model,levels = c("Observed", "ITT", "CACE")),
                           dv = factor(dv, levels = unique(model_effects$dv), labels = labels))

p <- ggplot(model_effects, aes(model, estimate)) + 
  geom_pointrange(aes(ymin = lo95, ymax = hi95), shape = 19, fatten = 1, size = 1) + 
  geom_pointrange(aes(ymin = lo83, ymax = hi83), shape = 19, fatten = 1, size = 2) + 
  geom_hline(yintercept = 0) + 
  facet_wrap( ~ dv, scales = "free", ncol = 3) +
  theme_ipsum_rc() + 
  theme(panel.background = element_blank(), 
        axis.title.y = element_blank(), 
        axis.title.x = element_blank(), 
        axis.text = element_text(size = 12), 
        plot.caption = element_text(size = 8), 
        strip.text.x = element_text(size = 12, colour = "darkblue", face = "bold"),
        plot.margin=unit(c(0.1,0.1,0.1,0.1), "cm"),
        panel.spacing = unit(0.25, "cm"))
p
# save plot
ggsave(p, file = "../output/wahlomat-other-effects.png", 
       width = 3.5, height = 4.5, units = "cm", dpi = 300, scale = 5)



## Combine knowledge model results; plot -------

model_effects <- rbind(model_effects_wahlomat_civic_knowledge,
                       model_effects_wahlomat_candidate_knowledge,
                       model_effects_wahlomat_event_knowledge
)

dat <- filter(model_effects, model == "CACE")
labels = paste(as.character(dat$dv), dat$scale)

model_effects <- transform(model_effects, 
                           model = factor(model,levels = c("Observed", "ITT", "CACE")),
                           dv = factor(dv, levels = unique(model_effects$dv), labels = labels))

p <- ggplot(model_effects, aes(model, estimate)) + 
  geom_pointrange(aes(ymin = lo95, ymax = hi95), shape = 19, fatten = 1, size = 1) + 
  geom_pointrange(aes(ymin = lo83, ymax = hi83), shape = 19, fatten = 1, size = 2) + 
  geom_hline(yintercept = 0) + 
  facet_wrap( ~ dv, scales = "free", nrow = 1) +
  theme_ipsum_rc() + 
  theme(panel.background = element_blank(), 
        axis.title.y = element_blank(), 
        axis.title.x = element_blank(), 
        axis.text = element_text(size = 12), 
        plot.caption = element_text(size = 8), 
        strip.text.x = element_text(size = 12, colour = "darkblue", face = "bold"),
        plot.margin=unit(c(0.1,0.1,0.1,0.1), "cm"),
        panel.spacing = unit(0.25, "cm"))
p
# save plot
ggsave(p, file = "../output/wahlomat-knowledge-effects.png",
       width = 3.5, height = 1.25, units = "cm", dpi = 300, scale = 5)





## MAIN: Knowledge on party positions -----------------------

dv <- "wahlomat_perceivedpos_correct_avg"
label <- "Knowledge on parties'\nissue positions"
hist(waves_df_wide[,dv])
scale <- "[0-1]"
model_effects_wahlomat_perceivedpos_correct_avg <- naive_itt_cace_estimation(dv, label, scale)

## BY ITEM: Knowledge on party positions -----------------------
# DIESEL
waves_df_wide$wahlomat_diesel_correct <- rowMeans(waves_df_wide[,c(
  "wahlomat_diesel_cdsu_correct", "wahlomat_diesel_spd_correct",
  "wahlomat_diesel_lin_correct", "wahlomat_diesel_gru_correct",
  "wahlomat_diesel_fdp_correct", "wahlomat_diesel_afd_correct"
)], na.rm=TRUE)
hist(waves_df_wide$wahlomat_diesel_correct)
summary(waves_df_wide$wahlomat_diesel_correct)

dv <- "wahlomat_diesel_correct"
label <- "Diesel"
scale <- "[0-1]"
model_effects_wahlomat_diesel <- naive_itt_cace_estimation(dv, label, scale)

waves_df_wide$wahlomat_fakenews_correct <- rowMeans(waves_df_wide[,c(
  "wahlomat_fakenews_cdsu_correct", "wahlomat_fakenews_spd_correct",
  "wahlomat_fakenews_lin_correct", "wahlomat_fakenews_gru_correct",
  "wahlomat_fakenews_fdp_correct", "wahlomat_fakenews_afd_correct"
)], na.rm=TRUE)
hist(waves_df_wide$wahlomat_fakenews_correct)
summary(waves_df_wide$wahlomat_fakenews_correct)

dv <- "wahlomat_fakenews_correct"
label <- "Fake news"
scale <- "[0-1]"
model_effects_wahlomat_fakenews <- naive_itt_cace_estimation(dv, label, scale)

waves_df_wide$wahlomat_healthcare_correct <- rowMeans(waves_df_wide[,c(
  "wahlomat_healthcare_cdsu_correct", "wahlomat_healthcare_spd_correct",
  "wahlomat_healthcare_lin_correct", "wahlomat_healthcare_gru_correct",
  "wahlomat_healthcare_fdp_correct", "wahlomat_healthcare_afd_correct"
)], na.rm=TRUE)
hist(waves_df_wide$wahlomat_healthcare_correct)
summary(waves_df_wide$wahlomat_healthcare_correct)

dv <- "wahlomat_healthcare_correct"
label <- "Health care"
scale <- "[0-1]"
model_effects_wahlomat_healthcare <- naive_itt_cace_estimation(dv, label, scale)

waves_df_wide$wahlomat_housing_correct <- rowMeans(waves_df_wide[,c(
  "wahlomat_housing_cdsu_correct", "wahlomat_housing_spd_correct",
  "wahlomat_housing_lin_correct", "wahlomat_housing_gru_correct",
  "wahlomat_housing_fdp_correct", "wahlomat_housing_afd_correct"
)], na.rm=TRUE)
hist(waves_df_wide$wahlomat_housing_correct)
summary(waves_df_wide$wahlomat_housing_correct)

dv <- "wahlomat_housing_correct"
label <- "Housing"
scale <- "[0-1]"
model_effects_wahlomat_housing <- naive_itt_cace_estimation(dv, label, scale)

waves_df_wide$wahlomat_renewable_correct <- rowMeans(waves_df_wide[,c(
  "wahlomat_renewable_cdsu_correct", "wahlomat_renewable_spd_correct",
  "wahlomat_renewable_lin_correct", "wahlomat_renewable_gru_correct",
  "wahlomat_renewable_fdp_correct", "wahlomat_renewable_afd_correct"
)], na.rm=TRUE)
hist(waves_df_wide$wahlomat_renewable_correct)
summary(waves_df_wide$wahlomat_renewable_correct)

dv <- "wahlomat_renewable_correct"
label <- "Renewable"
scale <- "[0-1]"
model_effects_wahlomat_renewable <- naive_itt_cace_estimation(dv, label, scale)

waves_df_wide$wahlomat_asylumlimit_correct <- rowMeans(waves_df_wide[,c(
  "wahlomat_asylumlimit_cdsu_correct", "wahlomat_asylumlimit_spd_correct",
  "wahlomat_asylumlimit_lin_correct", "wahlomat_asylumlimit_gru_correct",
  "wahlomat_asylumlimit_fdp_correct", "wahlomat_asylumlimit_afd_correct"
)], na.rm=TRUE)
hist(waves_df_wide$wahlomat_asylumlimit_correct)
summary(waves_df_wide$wahlomat_asylumlimit_correct)

dv <- "wahlomat_asylumlimit_correct"
label <- "Asylum limit"
scale <- "[0-1]"
model_effects_wahlomat_asylumlimit <- naive_itt_cace_estimation(dv, label, scale)


## Combine model results; plot -------

model_effects <- rbind(model_effects_wahlomat_diesel,
                       model_effects_wahlomat_fakenews,
                       model_effects_wahlomat_healthcare,
                       model_effects_wahlomat_housing,
                       model_effects_wahlomat_renewable,
                       model_effects_wahlomat_asylumlimit
)

dat <- filter(model_effects, model == "CACE")
labels = paste(as.character(dat$dv), dat$scale)

model_effects <- transform(model_effects, 
                           model = factor(model,levels = c("Observed", "ITT", "CACE")),
                           dv = factor(dv, levels = unique(model_effects$dv), labels = labels))

p <- ggplot(model_effects, aes(model, estimate)) + 
  geom_pointrange(aes(ymin = lo95, ymax = hi95), shape = 19, fatten = 1, size = 1) + 
  geom_pointrange(aes(ymin = lo83, ymax = hi83), shape = 19, fatten = 1, size = 2) + 
  geom_hline(yintercept = 0) + 
  facet_wrap( ~ dv, scales = "free", nrow = 2) +
  theme_ipsum_rc() + 
  theme(panel.background = element_blank(), 
        axis.title.y = element_blank(), 
        axis.title.x = element_blank(), 
        axis.text = element_text(size = 12), 
        plot.caption = element_text(size = 8), 
        strip.text.x = element_text(size = 12, colour = "darkblue", face = "bold"),
        plot.margin=unit(c(0.1,0.1,0.1,0.1), "cm"),
        panel.spacing = unit(0.25, "cm"))
p
# save plot
ggsave(p, file = "../output/wahlomat-knowledge-effects-by-item.png",
       width = 3.5, height = 2.5, units = "cm", dpi = 300, scale = 5)



